Research Question:

What is the effect of the teacher intervention on student grades among latinx students relative to their white peers?

Dataset structural features & theoretical considerations

Model levels:

  1. Grade-level (i.e., multiple course grades are nested within student clusters)
  2. Student-level
  3. Teacher-level

LEVEL REMOVED: Data has been averaged over courses within teacher clusters (i.e., one row per student-by-teacher-by-cohort). Courses is considered a nuisance-level in this research context especially because course subject data is not available (e.g., science, math, English, etc.).

We tried the IV modeling approach, but this method relies on a strong correlation between assignment and treatment. When the correlation between assignment and treatment is small, the instrument is said to be weak, which results in an invalid (biased) estimate of the causal effect. Therefore, other modeling approaches were considered.

2019 COHORT DROPPED: This cohort was dropped based on the significant environmental events effecting grades in 2019. It was deemed that the grading policy changes resulted in grades with qualitatively distinct meaning relative to previous cohorts making grade comparisons across cohorts unfeasible.

FOCAL VARIABLES:

CONTROL VARIABLES (covariates):

STUDENT-LEVEL

TEACHER-LEVEL

PRIOR (REMOVED): Prior year average grades have been removed from analyses as a covariate. THis is due to challenges with access to this data as well as concerns with the prior year scores being contaminated by treatment.

————————————————————————————–

## NULL

View all levels of course_subject

tabyl(droplevels(full_sample$course_subject))

View dataframe

library(DT)
datatable(full_sample, rownames = FALSE, filter="top",
          options = list(pageLength = 15, scrollX=T) )

Look at Treatment dose at student-level (prop_treat)

##  full_sample$prop_treat    n      percent
##               0.0000000 4537 0.2214142794
##               0.1666667   42 0.0020496803
##               0.2000000  285 0.0139085452
##               0.2500000 1036 0.0505587819
##               0.2857143   28 0.0013664536
##               0.3333333 2232 0.1089258699
##               0.3750000   16 0.0007808306
##               0.4000000  640 0.0312332243
##               0.4285714   28 0.0013664536
##               0.5000000 4328 0.2112146796
##               0.5555556    9 0.0004392172
##               0.5714286   91 0.0044409741
##               0.6000000  725 0.0353813870
##               0.6250000   40 0.0019520765
##               0.6666667 1905 0.0929676443
##               0.7142857  168 0.0081987214
##               0.7500000  716 0.0349421697
##               0.8000000  425 0.0207408130
##               0.8333333  228 0.0111268362
##               0.8571429   77 0.0037577473
##               0.8750000   24 0.0011712459
##               1.0000000 2911 0.1420623688

Create all group samples

See distribution of the 3-category course type variable

tabyl(full_sample$course_type3)
##  full_sample$course_type3    n   percent
##                  non_acad 5212 0.2543556
##                  non_stem 7873 0.3842175
##                      stem 7406 0.3614270

————————————————————————————–

Look at covariate balance

————————————————————————————–

Variable summary grouped by treatment status

full_sample %>%
   dplyr::select(treat, cohort, mark1, course_type3,
                 teach_age, teach_white_latinx, teach_gender, school_type, number_students_16, 
                 white_latinx, title, student_age, total_service_avg, sed, avg_absences) %>%
   tbl_summary(by = treat,
    statistic = list(all_continuous() ~ "{mean} ({sd})",
                     all_categorical() ~ "{n} ({p}%)")) %>%
  add_n() %>% add_p() %>%
  modify_header(label = "**Variable**") %>%
  bold_labels()
Variable N Control, N = 110061 Treatment, N = 94851 p-value2
cohort 20491 <0.001
2017 7682 (70%) 7709 (81%)
2018 3324 (30%) 1776 (19%)
mark1 20491 8.2 (3.4) 8.5 (3.3) <0.001
course_type3 20491 <0.001
non_acad 2167 (20%) 3045 (32%)
non_stem 3426 (31%) 4447 (47%)
stem 5413 (49%) 1993 (21%)
teach_age 20491 39 (11) 43 (11) <0.001
teach_white_latinx 20008 <0.001
White 7579 (71%) 7370 (79%)
Latinx 1118 (10%) 829 (8.9%)
Other 2013 (19%) 1099 (12%)
Unknown 296 187
teach_gender 19993 0.066
Male 3568 (33%) 3208 (35%)
Female 7142 (67%) 6075 (65%)
Unknown 296 202
school_type 20491 <0.001
JH 3881 (35%) 3893 (41%)
HS 7125 (65%) 5592 (59%)
number_students_16 17039 26 (16) 21 (11) <0.001
Unknown 2354 1098
white_latinx 18709 <0.001
White 4160 (41%) 3971 (46%)
Latinx 5920 (59%) 4658 (54%)
Unknown 926 856
title 20491 <0.001
Alt. Teacher 517 (4.7%) 206 (2.2%)
Teacher 10489 (95%) 9279 (98%)
student_age 20489 13.92 (1.78) 14.00 (1.82) 0.008
Unknown 0 2
total_service_avg 20491 12 (9) 16 (9) <0.001
sed 20489 <0.001
Not SED 5865 (53%) 5508 (58%)
SED 5141 (47%) 3975 (42%)
Unknown 0 2
avg_absences 20489 4.1 (4.5) 4.0 (4.5) 0.2
Unknown 0 2

1 Statistics presented: n (%); mean (SD)

2 Statistical tests performed: chi-square test of independence; Wilcoxon rank-sum test

————————————————————————————–

Multi-level modeling

————————————————————————————–

library(psych)
library(rockchalk)
library(lme4)
library(varTestnlme)
library(lmerTest)
library(clubSandwich)
library(effects)

options(scipen = 999)

source(here("R-MLM-Functions.R"))

————————————————————————————–

STEM COURSE SAMPLE (Teacher-level N = 65)

  • treat: has variation at student & teacher level

  • white_latinx (level2): center within teacher?

    • varies at the teacher-level (proportion of Latinx varies within teacher cluster)
  • If white_latinx is group-mean centered at the teacher-level could allow the between-level variable, within-level variable, or both to interact with treatment.

————————————————————————————–

Centering

  • cwc: group-mean centered (pure within)
  • mj_cgm: grand-mean centered (pure between)
  • cgm: grand-mean centered (smushed)
stem_prepared <- stem_sample %>% 
  mutate(white_latinx = as.numeric(white_latinx) - 1) %>% 
  mutate(cohort = as.numeric(cohort) - 1) %>% 
  mutate(school_type = as.numeric(school_type) - 1)


stem_centered <- gmc(stem_prepared,
                 c("white_latinx"), 
                 by = c("teacher_id"),
                 FUN = mean,
                 suffix = c("_meanj", "_cwc"),
                 fulldataframe = TRUE) %>% 
mutate(white_latinx_mj_cgm = white_latinx_meanj - mean(white_latinx_meanj, na.rm = TRUE)) %>% 
mutate(cohort_cgm = cohort - mean(cohort, na.rm = TRUE)) %>% 
mutate(school_type_cgm = school_type - mean(school_type, na.rm = TRUE))

————————————————————————————–

Empty model with teacher-level & student-level random intercepts

mlm_level23 <- lmer(mark1 ~ 1 + (1 | teacher_id) + (1 | student_id), data = stem_sample, REML = FALSE)
mlm_level23 
## Linear mixed model fit by maximum likelihood  ['lmerModLmerTest']
## Formula: mark1 ~ 1 + (1 | teacher_id) + (1 | student_id)
##    Data: stem_sample
##       AIC       BIC    logLik  deviance  df.resid 
##  37144.11  37171.75 -18568.06  37136.11      7402 
## Random effects:
##  Groups     Name        Std.Dev.
##  student_id (Intercept) 2.313   
##  teacher_id (Intercept) 1.568   
##  Residual               2.050   
## Number of obs: 7406, groups:  student_id, 5330; teacher_id, 65
## Fixed Effects:
## (Intercept)  
##       7.225
# ICC = (VAR_B1 + VAR_B2) / (VAR_B1 + VAR_B2 + VAR_W)
mlm <- as.data.frame(VarCorr(mlm_level23))

icc_students <-  ((mlm[1,4] + mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))

icc_teachers <-  ((mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))

# Proportion of variance captured by `student-level` random intercept:
round(icc_students*100,2)
## [1] 64.99
# Proportion of variance captured by `teacher-level` random intercept:
round(icc_teachers*100,2)
## [1] 20.46

Uncentered white_latinx model

mlm_stem <- lmer(mark1 ~ 1 +
                   treat*white_latinx +
                   cohort_cgm + 
                   school_type_cgm +
               (1 | teacher_id) + (1 | student_id), data = stem_centered)

summary(mlm_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx + cohort_cgm + school_type_cgm +  
##     (1 | teacher_id) + (1 | student_id)
##    Data: stem_centered
## 
## REML criterion at convergence: 33440.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03092 -0.46349  0.06627  0.50531  2.77731 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  student_id (Intercept) 4.942    2.223   
##  teacher_id (Intercept) 1.358    1.165   
##  Residual               4.059    2.015   
## Number of obs: 6744, groups:  student_id, 4874; teacher_id, 65
## 
## Fixed effects:
##                              Estimate Std. Error        df t value
## (Intercept)                    8.5516     0.1909   67.4460  44.796
## treatTreatment                -0.4653     0.3503   65.7620  -1.328
## white_latinx                  -2.2317     0.0973 5752.1595 -22.937
## cohort_cgm                    -0.4566     0.3284   56.5169  -1.390
## school_type_cgm               -1.2682     0.3382   56.1374  -3.749
## treatTreatment:white_latinx    0.5954     0.1635 4784.5016   3.641
##                                         Pr(>|t|)    
## (Intercept)                 < 0.0000000000000002 ***
## treatTreatment                          0.188694    
## white_latinx                < 0.0000000000000002 ***
## cohort_cgm                              0.169841    
## school_type_cgm                         0.000421 ***
## treatTreatment:white_latinx             0.000274 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtTrt wht_lt chrt_c schl__
## treatTrtmnt -0.526                            
## white_latnx -0.296  0.110                     
## cohort_cgm   0.015 -0.105 -0.010              
## schl_typ_cg  0.020 -0.106  0.009 -0.203       
## trtTrtmnt:_  0.123 -0.279 -0.417 -0.011 -0.026

Centered white_latinx within & between components interact

mlm_stem <- lmer(mark1 ~ 1 +
                   treat*white_latinx_cwc + 
                   treat*white_latinx_mj_cgm + 
                   cohort_cgm + 
                   school_type_cgm +
               (1 | teacher_id) + (1 | student_id), data = stem_centered)

summary(mlm_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx_cwc + treat * white_latinx_mj_cgm +  
##     cohort_cgm + school_type_cgm + (1 | teacher_id) + (1 | student_id)
##    Data: stem_centered
## 
## REML criterion at convergence: 33411.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.94604 -0.46685  0.06255  0.50382  2.76688 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  student_id (Intercept) 4.8740   2.2077  
##  teacher_id (Intercept) 0.8274   0.9096  
##  Residual               4.1050   2.0261  
## Number of obs: 6744, groups:  student_id, 4874; teacher_id, 65
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                           7.49198    0.14904   51.85344  50.267
## treatTreatment                       -0.08519    0.27628   49.82792  -0.308
## white_latinx_cwc                     -2.19194    0.09762 5717.32442 -22.454
## white_latinx_mj_cgm                  -5.50521    0.72585   54.22499  -7.585
## cohort_cgm                           -0.21044    0.26727   53.02994  -0.787
## school_type_cgm                      -1.25513    0.27743   52.79610  -4.524
## treatTreatment:white_latinx_cwc       0.62027    0.16494 4772.85664   3.761
## treatTreatment:white_latinx_mj_cgm    0.56014    1.17533   54.50462   0.477
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## treatTreatment                                 0.759095    
## white_latinx_cwc                   < 0.0000000000000002 ***
## white_latinx_mj_cgm                      0.000000000452 ***
## cohort_cgm                                     0.434589    
## school_type_cgm                          0.000034739322 ***
## treatTreatment:white_latinx_cwc                0.000171 ***
## treatTreatment:white_latinx_mj_cgm             0.635569    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtTrt wht_l_ wht___ chrt_c schl__ trT:__
## treatTrtmnt -0.529                                          
## wht_ltnx_cw -0.002 -0.003                                   
## wht_ltnx_m_ -0.203  0.107  0.017                            
## cohort_cgm   0.041 -0.096  0.003 -0.124                     
## schl_typ_cg -0.001 -0.058 -0.002  0.138 -0.193              
## trtTrtmn:__  0.001  0.001 -0.416 -0.003 -0.001 -0.001       
## trtTrtm:___  0.116 -0.222  0.002 -0.615 -0.031 -0.223  0.002

————————————————————————————–

NON-STEM COURSE SAMPLE (Teacher-level N = 85)

————————————————————————————–

non_stem_prepared <- non_stem_sample %>% 
  mutate(white_latinx = as.numeric(white_latinx) - 1) %>% 
  mutate(cohort = as.numeric(cohort) - 1) %>% 
  mutate(school_type = as.numeric(school_type) - 1)

non_stem_centered <- gmc(non_stem_prepared,
                 c("white_latinx"), 
                 by = c("teacher_id"),
                 FUN = mean,
                 suffix = c("_meanj", "_cwc"),
                 fulldataframe = TRUE) %>% 
mutate(white_latinx_mj_cgm = white_latinx_meanj - mean(white_latinx_meanj, na.rm = TRUE)) %>% 
mutate(cohort_cgm = cohort - mean(cohort, na.rm = TRUE)) %>% 
mutate(school_type_cgm = school_type - mean(school_type, na.rm = TRUE))

Empty model with teacher-level & student-level random intercepts

mlm_level23 <- lmer(mark1 ~ 1 + (1 | teacher_id) + (1 | student_id), data = non_stem_sample, REML = FALSE)
mlm_level23 
## Linear mixed model fit by maximum likelihood  ['lmerModLmerTest']
## Formula: mark1 ~ 1 + (1 | teacher_id) + (1 | student_id)
##    Data: non_stem_sample
##       AIC       BIC    logLik  deviance  df.resid 
##  38666.95  38694.83 -19329.47  38658.95      7869 
## Random effects:
##  Groups     Name        Std.Dev.
##  student_id (Intercept) 2.416   
##  teacher_id (Intercept) 1.436   
##  Residual               1.772   
## Number of obs: 7873, groups:  student_id, 5775; teacher_id, 80
## Fixed Effects:
## (Intercept)  
##       7.917
# ICC = (VAR_B1 + VAR_B2) / (VAR_B1 + VAR_B2 + VAR_W)
mlm <- as.data.frame(VarCorr(mlm_level23))

icc_students <-  ((mlm[1,4] + mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))

icc_teachers <-  ((mlm[2,4]) / (mlm[1,4] + mlm[2,4] + mlm[3,4]))

# Proportion of variance captured by `student-level` random intercept:
round(icc_students*100,2)
## [1] 71.56
# Proportion of variance captured by `teacher-level` random intercept:
round(icc_teachers*100,2)
## [1] 18.67

Uncentered white_latinx model

mlm_non_stem <- lmer(mark1 ~ 1 +
                   treat*white_latinx +
                   cohort_cgm + 
                   school_type_cgm +
               (1 | teacher_id) + (1 | student_id), data = non_stem_centered)

summary(mlm_non_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx + cohort_cgm + school_type_cgm +  
##     (1 | teacher_id) + (1 | student_id)
##    Data: non_stem_centered
## 
## REML criterion at convergence: 34865.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8460 -0.3878  0.0840  0.4588  3.8396 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  student_id (Intercept) 5.235    2.288   
##  teacher_id (Intercept) 1.622    1.273   
##  Residual               3.093    1.759   
## Number of obs: 7188, groups:  student_id, 5300; teacher_id, 80
## 
## Fixed effects:
##                               Estimate Std. Error         df t value
## (Intercept)                    9.17439    0.23919   94.31053  38.356
## treatTreatment                -0.10874    0.31938   86.66579  -0.340
## white_latinx                  -2.01075    0.11646 7097.65153 -17.265
## cohort_cgm                     0.01963    0.33671   73.58712   0.058
## school_type_cgm               -0.97013    0.32118   75.24210  -3.021
## treatTreatment:white_latinx   -0.07065    0.13851 6083.05070  -0.510
##                                         Pr(>|t|)    
## (Intercept)                 < 0.0000000000000002 ***
## treatTreatment                           0.73432    
## white_latinx                < 0.0000000000000002 ***
## cohort_cgm                               0.95368    
## school_type_cgm                          0.00345 ** 
## treatTreatment:white_latinx              0.61001    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtTrt wht_lt chrt_c schl__
## treatTrtmnt -0.748                            
## white_latnx -0.312  0.195                     
## cohort_cgm  -0.133  0.151 -0.004              
## schl_typ_cg -0.168  0.174  0.004 -0.062       
## trtTrtmnt:_  0.223 -0.267 -0.704 -0.015  0.000

Centered white_latinx (interaction for both within & between components)

mlm_non_stem <- lmer(mark1 ~ 1 +
                   treat*white_latinx_cwc + 
                   treat*white_latinx_mj_cgm + 
                   cohort_cgm + 
                   school_type_cgm +
               (1 | teacher_id) + (1 | student_id), data = non_stem_centered)

summary(mlm_non_stem)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mark1 ~ 1 + treat * white_latinx_cwc + treat * white_latinx_mj_cgm +  
##     cohort_cgm + school_type_cgm + (1 | teacher_id) + (1 | student_id)
##    Data: non_stem_centered
## 
## REML criterion at convergence: 34853.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8338 -0.3886  0.0828  0.4595  3.8011 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  student_id (Intercept) 5.221    2.285   
##  teacher_id (Intercept) 1.476    1.215   
##  Residual               3.101    1.761   
## Number of obs: 7188, groups:  student_id, 5300; teacher_id, 80
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                           8.00568    0.22821   71.85462  35.081
## treatTreatment                       -0.13977    0.30274   70.89531  -0.462
## white_latinx_cwc                     -2.00775    0.11696 7044.93520 -17.166
## white_latinx_mj_cgm                  -2.24742    0.99173   75.95655  -2.266
## cohort_cgm                            0.29689    0.33649   71.97193   0.882
## school_type_cgm                      -1.01251    0.30798   73.03861  -3.288
## treatTreatment:white_latinx_cwc      -0.04858    0.13918 6013.27890  -0.349
## treatTreatment:white_latinx_mj_cgm   -2.69726    1.40016   74.70703  -1.926
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## treatTreatment                                  0.64573    
## white_latinx_cwc                   < 0.0000000000000002 ***
## white_latinx_mj_cgm                             0.02629 *  
## cohort_cgm                                      0.38054    
## school_type_cgm                                 0.00156 ** 
## treatTreatment:white_latinx_cwc                 0.72709    
## treatTreatment:white_latinx_mj_cgm              0.05786 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtTrt wht_l_ wht___ chrt_c schl__ trT:__
## treatTrtmnt -0.764                                          
## wht_ltnx_cw -0.003 -0.002                                   
## wht_ltnx_m_ -0.296  0.225  0.020                            
## cohort_cgm  -0.153  0.162 -0.008  0.046                     
## schl_typ_cg -0.172  0.180  0.002  0.022 -0.071              
## trtTrtmn:__  0.002  0.002 -0.705 -0.006  0.007 -0.001       
## trtTrtm:___  0.234 -0.182 -0.002 -0.714 -0.232  0.014  0.000

Show NON-STEM & STEM model result side-by-side

library(sjPlot)  # Latex MLM tables

tab_model(mlm_non_stem,mlm_stem,
          dv.labels = c("NON-STEM","STEM"),
          show.se = TRUE,
          show.ci = FALSE
          #p.val = "kr",    # Kenward-Roger correction
          #show.df = TRUE
          )
  NON-STEM STEM
Predictors Estimates std. Error p Estimates std. Error p
(Intercept) 8.01 0.23 <0.001 7.49 0.15 <0.001
treat [Treatment] -0.14 0.30 0.644 -0.09 0.28 0.758
white_latinx_cwc -2.01 0.12 <0.001 -2.19 0.10 <0.001
white_latinx_mj_cgm -2.25 0.99 0.023 -5.51 0.73 <0.001
cohort_cgm 0.30 0.34 0.378 -0.21 0.27 0.431
school_type_cgm -1.01 0.31 0.001 -1.26 0.28 <0.001
treat [Treatment] *
white_latinx_cwc
-0.05 0.14 0.727 0.62 0.16 <0.001
treat [Treatment] *
white_latinx_mj_cgm
-2.70 1.40 0.054 0.56 1.18 0.634
Random Effects
σ2 3.10 4.11
τ00 5.22 student_id 4.87 student_id
1.48 teacher_id 0.83 teacher_id
ICC 0.68 0.58
N 80 teacher_id 65 teacher_id
5300 student_id 4874 student_id
Observations 7188 6744
Marginal R2 / Conditional R2 0.144 / 0.729 0.202 / 0.666